Reptational dynamics in dissipative particle dynamics simulations of polymer melts 



Petri Nikunen 

Biophysics and Statistical Mechanics Group, Laboratory of Computational Engineering, 
Helsinki University of Technology, P.O. Box 9203, Fl-02015 HUT, Finland 



in 
o 
o 

(N 
O 

Q 

(N 



o 

CO 

-4-* 



Ilpo Vattulainen 

Laboratory of Physics and Helsinki Institute of Physics, Helsinki University of Technology, P.O. Box 1100, 
Fl-02015 HUT, Finland; M EM PHYS-C 'enter for Biomembrane Physics, University of Southern Denmark, 
Odense, Denmark; Institute of Physics, Tampere University of Technology, P.O. Box 692, FI-33101 Tampere, Finland 

Mikko Karttunen 

Biophysics and Statistical Mechanics Group, Laboratory of Computational Engineering, 
Helsinki University of Technology, P.O. Box 9203, Fl-02015 HUT, Finland; 
Department of Applied Mathematics, University of Western Ontario, London ( ON), Canada 

(Dated: November 27, 2005) 

Understanding the complex viscoelastic properties of polymeric liquids remains a challenge in materials sci- 
ence and soft matter physics. Here, we present a simple and computationally efficient criterion for the topolog- 
ical constraints in polymeric liquids using the Dissipative Particle Dynamics (DPD). The same approach is also 
applicable in other soft potential models. For short chains the model correctly reproduces Rouse-like dynamics 
whereas for longer chains the dynamics becomes reptational as the chain length is increased — something that is 
not attainable using standard DPD or other coarse-grained soft potential methods. Importantly, no new length 
scales or forces need to be added. 
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I. INTRODUCTION 

The static and dynamic properties of polymeric liquids are, 
by and large, dominated by topological constraints. The ori- 
gin of these constraints is easy to understand: polymers can 
slide past but not penetrate through each other. That is the 
physical origin of the reptation model 0, IH which is the 
most successful theory in describing the behavior of entan- 
gled polymers. Despite active research in the field, entan- 
gled polymeric liquids keep posing many challenges to the- 
orists U EH EL e xperimentalists fzl IH, @] and computational 
modelers 1 HJ, HH LL2, LLJ, LLJ, LL2l 1 1 611 . The importance of un- 
derstanding the fundamentals of polymeric liquids can hardly 
be overemphasized as they are one of the key issues in novel 
(bio)materials science I171M 811 . 

The dynamics of polymer melts is typically described in 
terms of the Rouse and reptation models |3]. Short chains are 
able to move to any direction and are not entangled. That is the 
physical origin of the Rouse model . For longer chains, 

entanglements and uncrossability of chains cannot be ignored, 
and the chains become constrained to move in the direction of 
the chain backbone. This motion resembles that of a reptating 
snake — hence the name reptation model llj|2j,|3j]. 

Computer simulations offer a detailed look into polymers 
and their dynamics. In classical molecular dynamics simula- 
tions the system size and simulation time pose limits as they 
are typically of the order of lOnm in linear size and around 
10 ns in time. In contrast, coarse grained methods, such as 
dissipative particle dynamics (DPD), allow access to microm- 
eter and microsecond scales. That is due to the soft potentials, 
and, like everything in life, they do not come without a price 
to pay: the softness of the conservative potentials allows the 
chains to slide through each other thus strongly affecting the 
dynamics of the system. Indeed, the scaling laws obtained 



from DPD simulations of polymer melts 1 20, 21] are not able 
to describe entangled liquids. This is a direct consequence 
of the fact that in DPD simulations, polymers can penetrate 
through themselves. Whereas that is not a problem in study- 
ing the equilibrium properties in the Rouse regime, reptation 
cannot be studied using the basic DPD model with soft inter- 
actions. 

To preserve the advantages of coarse-grained models and to 
correct for their deficiencies, Padding and Briels 1 13] recently 
introduced an algorithm that explicitly detects and prevents 
bond crossings. They consider bonds as elastic bands that be- 
come entangled and use energy minimization to determine the 
entanglement positions. This approach is general and very 
promising but it is also complicated to implement and compu- 
tationally intensive fkUl . 

Another promising approach was put forward by Pan and 
Manke |22]. They reduce the frequency of bond crossings by 
introducing segmental repulsive forces between the points of 
nearest contact between neighboring chains. This approach is 
simple to implement but the introduction of a new force and 
a related cutoff increases the computational load, and adds a 
new length scale whose physical determination is somewhat 
ambiguous. On the other hand, the model seems to be able to 
capture both the Rouse and reptational behavior l22ll like that 
of Padding and Briels ifllEl . 

In this article, we introduce a simple and generic criterion 
based on simple geometrical arguments to solve the crossabil- 
ity problem. No new forces are added, the approach is con- 
ceptually simple and does not depend on the level of coarse 
graining. Importantly, it allows easy, and if necessary even 
on-the-fiy, tuning between the Rouse and reptation regimes. 

The rest of this paper is organized as follows. In the next 
section we will briefly describe the DPD method. Section^] 
describes our criterion for including topological constraints in 
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DPD, or for that matter any other soft potential, simulation. 
In Sec.[V]\ve show results from our simulations and compare 
them to other methods. Finally, we finish with a discussion 
and outlook in Sec. lVll 



II. DISSIPATIVE PARTICLE DYNAMICS 

In DPD, the time evolution of particles is given by the New- 
ton's equations of motion, and the total force acting on particle 
i is given as a sum of pairwise conservative, dissipative, and 
random forces, respectively, as F t — (-^y + ^y + ^y ) ■ 

The conservative force is independent of the dissipative and 
random forces. Typically it takes the form 



ay(l - nj/r c )eij, if ry < r c 
0, otherwise, 



(1) 



with r 



— rj, ry = |ry|, andey = fij/rij. The variable 
describes the repulsion between particles i and j, and thus 



produces excluded volume interactions. 
The dissipative force is expressed as 



F. 



(2) 



where 7 is a friction parameter, uj D (rij) a weight function 



for the dissipative force, and Vi 



The dissipative 



force slows down the particles by decreasing kinetic energy 
from them. This effect is balanced by the random force due to 
thermal fluctuations, 




FIG. 1 : Two bonds crossing each other. If Eq. J5J is satisfied, cross- 
ings cannot occur. 



III. TOPOLOGICAL CONSTRAINTS 



To take into account the topological constraints, chain 
crossings must be prevented. As discussed in the introduc- 
tion, there are currently two off-lattice methods fl3l 12211 for 
this purpose. Here, we introduce a third alternative. 

First, each individual bead has a radius r m ; n which is im- 
penetrable to other beads. In systems with Lennard-Jones 
potentials that condition is automatically satisfied due to the 
r~ 12 part that takes care of the Fermi exclusion principle. 
In mesoscopic simulations with soft potentials that constraint 
needs special attention. Second, the intramolecular bonds 
have some maximum stretch, ^ max . By using simple geom- 
etry, we can postulate that if the condition 



F, 



(TUJ R (rij% 



(3) 



V2r 



> 



(5) 



where a is the amplitude of thermal noise, uj R (rij) is the 
weight function for the random force, and dj(t) are Gaus- 
sian random variables with (Cy(i)) = and (Cij{t)Ckl{t')) — 
(SikSji + Sii5 jk )S(t - if). The condition Cy(*) = Cji(*) is re- 
quired for momentum conservation. That is a necessary con- 
dition for the conservation of hydrodynamics. 

The weight functions u> D {rij) and u R (rn) cannot be cho- 
sen arbitrarily. Espanol and Warren 12411 showed that the 



fluctuation-dissipation relations u) (ry) 



- l.,R 



(ry)] 2 and 



a = 2jksT must be satisfied for the system to have a 
canonical equilibrium distribution. Here T is the temperature 
of the system and ku is the Boltzmann constant. The func- 
tional form of the weight functions is not defined by the DPD 
method but virtually all DPD studies use 



u D { nj ) = [cj R ( nj )} 2 = 



(1 - rij/r c ) 2 , if ry < r c 
0, otherwise. 



(4) 



Coarse graining in DPD comes in through the soft conser- 
vative potential and forces (Eq. Q). A detailed account of 
DPD, derivation of time and length scales, and its applications 
is given by R.D. Groot I25ll . An in depth discussion of coarse 
graining by P. Espanol can be found in the same reference. 



is satisfied, any two bonds cannot cross each other, see Fig.^ 
The length scales involved, i. e., r m i n and (t max have a clear 
physical meaning. 

The obvious question is whether that condition is actually 
useful and when does it work. As an example, let us con- 
sider DPD simulations of polymers. The parameters used in 
these simulations are often justified on the basis of the Flory- 
Huggins theory IzMErjll . where the key component is the solu- 
bility as expressed by the x-parameters. Then, in simulations 
of block co-polymers, e.g., it is the mutual repulsion between 
the different components that matters — as a matter of fact, 
the values of the interaction parameters ay may be derived 
in different ways and their values tell only about the degree 
of coarse-graining. The condition set by Eq. 10 can thus be 
met by a proper degree of coarse graining, complemented by 
a reasonable description for bond stretching (springs). Indeed, 
above £ max is limited by the type of springs used in the model. 
With FENE springs 1 27] that is easy to tune as they have only 
finite extension after which the force becomes infinite. With 
harmonic springs more care is needed to satisfy Eq. (|5ji as 
there is no FENE-like cutoff set by the equation of motion. 
We will return to that in the Sec.lVl 
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a (amplitude of the conservative force) 


50 


100 


150 


200 


k (spring constant) 


100 


200 


300 


400 


At (time step) 


0.03 


0.02 


0.015 


0.0125 



TABLE I: The parameters used in this study. 



IV. SIMULATIONS 

For simplicity, and to be able to compare the model with 
other simulations, we considered a melt of linear polymers 
in a cubic box (3D) with periodic boundary conditions. To 
avoid finite size effects, the linear box size L was chosen to be 
at least 1.75 times the average end-to-end distance of chains. 
We also carried out simulations with different box sizes to 
ensure that the systems were free of finite size effects. That 
was done since it is known that static properties are affected 
by them O . 

All the systems had 128 chains consisting of N monomers, 
and no additional solvent or free monomers were present. All 
monomers were chosen to be identical, and thus the monomer 
mass was set equal to unity, m = 1, fixing the scale of mass. 
The cutoff distance r c sets the length scale for the model. The 
conservative forces had the form given in Eq. Q, with r c = 1 
and dij = a for all particle pairs. The values of a, as well 
as other simulation parameters used in these simulations are 
listed in TablellVl 

For the random and dissipative forces we used Eqs. (|2ji and 
with the common choices 1251 12(1 l29t I3CM1 7 = 4.5 and 
(7 = 3. This sets the temperature to fc^T = 1, and hence the 
time scale is given by ^Jmr^Jk^T . 

The monomers were connected using harmonic springs, 
i.e., Ff = J2j k(£ — Tij)eij, where the sum runs over all 
particles j to which particle i is connected. The equilibrium 
bond length was set to £ = 0.95. That particular value was 
chosen as it is very near the first maximum of the radial dis- 
tribution function (at the density p — 1). The spring constant 
was chosen to be k = 2a. If k is much smaller, bonds are very 
flexible and Eq. (JSJl is not satisfied. On the other hand, if k is 
much larger than a, the time step At must be decreased from 
the value set by the choice of a thus decreasing the computa- 
tional efficiency. Another possibility would be to use FENE 
springs |28] since they have finite extension. 

The density was chosen to be p = 1, which is lower than 
the densities typically used in DPD simulations j2^,E3l- The 
reason for high densities is to give different repulsive interac- 
tions for different particle types. This works only if particles 
overlap each other considerably. In the present work, we don't 
need such interactions, and therefore the lower density is suffi- 
cient. In fact, the density of p = 1 sets the monomer-monomer 
coordination number near 12, which is a typical value for real 
liquids. 

All systems were started from random flight initial config- 
urations and they were equilibrated for 10 6 time steps. After 
the equilibration, we simulated systems at least for 10 7 time 
steps to compute the desired quantities. Equations of motion 
were integrated using the DPD-VV algorithm !29l l3lll32Tl . 
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FIG. 2: The snapshots present 10 superpositions of configurations for 
a chain of length TV = 256 taken at times 100 apart. Green: times up 
to 500 (in DPD time units), red: times from 500 to 1000. a) Rouse 
dynamics (a — 25, k = 50) and b) reptation (a = 100, k = 200). 



V. RESULTS 

Figure |2 shows snapshots of the chain motion during the 
simulation at different times and regimes. For clarity, the 
chain is projected onto two dimensions. It is immediately 
clear that the motions in Figs. |2^ and |2j) are qualitatively 
different. Figure |2ji shows Rouse-like motion in which the 
polymers are free to move in every direction, and Fig. |2j? rep- 
resents reptation confined into a tube. 



A. Radial distribution function 

We will now study the static properties to see the effect of 
Eq. As discussed, by tuning the chain stiffness it is pos- 
sible to move gradually from the Rouse regime to reptation. 
This should be reflected in both the radial distribution func- 
tion and the bond length distribution. 

The radial distribution function (RDF) g(r) describes the 
qualitative structure of a fluid. It is defined as g(r) = p{r)/p 
where p(r ) is the average density from a given particle at a dis- 
tance r. Figures [3^ and|3j> show the radial distribution func- 
tion g(r) and the bond length distribution fi{r) for different 
parameter sets for chains of length N = 32. The arrows in the 
figures indicate the values of r m - ln and £ max in Eq. (JSJl. As the 
figures show, Eq. (|3J is satisfied for larger values of a and k. 
The small non-zero values below r m i n are due to the softness 
of the interparticle DPD potentials. 

The above can be characterized by taking a look at the 
average bond lengths are their mean square deviations. For 
dij = 50 we measured (£) = 0.977 ± 0.092. As the strength 
of interaction is increased we obtain (£) = 0.968 ± 0.064 
for dij = 100, {£) = 0.966 ± 0.051 for ay = 150, and 
{£) = 0.965 ± 0.045 for ay = 200. The most important 
issue is the decrease of the mean square deviation as that re- 
stricts the amount of overlap between the monomers of differ- 
ent chains. Importantly, for FENE chains this can be directly 
controlled by using the above measurements and RDF as a 
guideline and setting the maximum extent of the chain to an 
appropriate value. 
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FIG. 3: a) Radial distribution function in the case of N = 32. The 
arrow shows the distance r m i n , and the inset the region at lengths 
shorter than r m i B . Compare with LJ models (r c — 1, 2.5). b) Bond 
length distribution (TV = 32). The arrow shows the location of l m&K , 
and the inset the region at values larger than £ max . 




FIG. 4: a) The radius of gyration as a function of chain length for 
different parameters, b) The end-to-end vector autocorrelation (o = 
100, k = 200) for chains of different length. 



A comparison of the radial distribution functions shows 
that the current approach allows for tuning between typical 
DPD results Il2lll26l 13211 and typical molecular dynamics sim- 
ulations using Lennard- Jones potentials [ 281] . As the bond 
strength is increased (a > 100, k > 200), g(r) becomes qual- 
itatively similar to that from a Lennard- Jones system. 



B. Static scaling 

Next, we studied the end-to-end distance R and the radius 
of gyration R g . The former is defined as R — \R\ = | rt — f$}\ 

rim) 2 , where r cm = 



.iy 9 ( 



and the latter as R 2 g 

Sili Ti. In a polymer melt, they are expected to scale 



j_ 

JV 

as (R) oc N 1 /' 2 and (R g ) oc N 1 / 2 in both Rouse and the 
reptation regime. Previous studies using soft potentials 12011 
and systems with more realistic hard potentials 1 28] exhibit 
scaling. In Fig.@k we plot the results for the radius of gyration 
for different parameter sets. It is clear from the figure that the 
system exhibits the proper scaling behavior independently of 
the interaction parameters as it should. 




FIG. 5: a) Scaling of the longest relaxation time r. There is a 
crossover from Rouse scaling (r oc TV 2 ) to reptation (r oc N [i ). 
b) Similarly, the proper scaling limits are reached for the diffusion 
coefficient D. 



C. Relaxation time 

One of the main practical obstacles in simulations of poly- 
meric solutions is the long stress relaxation time. The longest 
relaxation time, t, depends on the molecular weight and the 
reptation theory predicts it to scale as r oc TV 3 . That pre- 
diction assumes only one mechanism for relaxation, i. e., dif- 
fusion along the countour [1J]. The Rouse model predicts a 
distinctly different behavior with r oc /V 2 . 

To estimate the scaling behavior, we measured the end-to- 
end autocorrelation function. It is shown in Fig. for poly- 
mers of different length. Assuming exponential decay, i. e., 

(R(t) ■ R(0)) ~ exp(-i/r), 

we can extract the longest relaxation time r by fitting. Fig- 
ure^ shows that both scaling regimes are captured properly. 
Fig. [5^ illustrates one of the main results of this paper: the 
simple criterion summarized by Eq. allows an easy, phys- 
ical and computationally efficient tuning between the Rouse 
regime and reptation. 

The scaling exponents 2 and 3 for Rouse and reptation, in 
respective order, are the limiting laws. The exponents have 
been frequently debated in the literature. For example, for 
the same values of /V as used here, Padding and Briels l23ll 
found two scaling regimes in t, one with exponent 2.8 and 
the other one with exponent 3.5. The dependence t oc N 3A is 
experimentally observed for the longest relaxation time in the 
entangled regime 13311 . This discrepancy is often associated 
with fluctuations of the contour length of the primitive path. In 
a real situation, however, the tube has a characteristic lifetime 
and the length of the primitive path fluctuates since the Rouse 
modes continue in the direction along the primitive path. 

Determining the value of the exponent was not the main 
goal here, and thus we did not attempt to evaluate it in a high 
precision — we will focus on that in a future publication. The 
above simply to demonstrates that it is indeed possible to use 
the soft DPD model to describe entangled polymer liquids 
realistically without introducing additional length scales and 
forces. 
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D. Diffusion 

The motion of a polymer, or its segments, is described by 
the diffusion coefficient. Typically, one measures the center- 
of-mass diffusion coefficient for a polymer chain, i. e., 

D= lim i([r cm (£)-f* cm (0)] 2 ). 

t— >oo D 

The scaling of D with molecular weight has been studied 
intensely over the years, see e. g. the discussion in £1- 

The theory predicts two scaling limits, D oc TV -1 for the 
Rouse model and D oc 7V~ 2 for the pure reptation model (N 
is proportional to molecular weight). Figure [3J> shows that 
both scaling regimes are found. 

Considering the nature of the DPD model, it is remarkable 
that both regimes are recovered. In simulations using the plain 
DPD model without paying attention to the criterion given by 
Eq. 0, only Rouse scaling has been observed I20l 12111 even 
in the case of long polymers. 

As with the longest relaxation time, the exponents typ- 
ically reported are between the scaling limits. Pearson et 
al. 13411 measured D as a function of molecular weight M w 
in polyethylene and they found that the diffusion coefficient 
follows a power law D = 1.65Af~ 1 ' 98 cm 2 /s for the entire 
range from M w = 600 to = 12000 (g/mol). The simula- 
tions by Kremer et al. OH and Padding and Briels lUlEl 
confirmed this finding: the center of mass diffusion coefficient 
scales as D oc 7V~ 2 in melt. 

Padding and Briels [23] compared their results to different 
simulations and experiments, and found that in ethylene the 
crossover between Rouse-like and reptational dynamics takes 
place at molecular weight of 560 g/mol (corresponds to 40 
ethylenes). Because in Fig. |5j> the crossover takes place be- 
tween N = 40 and N = 60, we can picture each particle 
roughly as one ethylene unit. 

VI. DISCUSSION 

In this article we have presented a simple criterion for topo- 
logical constraints in coarse grained DPD simulations of poly- 



meric liquids. No new forces or length scales were added. 
We showed that this approach is able to reproduce the Rouse 
model at one limit and reptational dynamics at the other. Here, 
we validated and demonstrated this approach against other 
models and experimental results using linear homopolymers. 
This approach can also be used for systems of, e. g., block co- 
polymers with different interactions and monomer sizes, and 
shear simulations. In practice, one can always run a short test 
simulation, and use g(r) and the bond length distribution (as 
in Figs.|3t and b) to verify that the criterion set by Eq. (0 is 
met. 

There is one other issue that we need to address, namely 
by tuning the chain stiffness one inevitably changes the en- 
tanglement length in addition to intercrossability of chains. 
It is known from previous simulations using Lennard-Jones 
as well as some coarse-grained models that increasing chain 
stiffness intensifies reptation I 35ll36[l37ll . To account for this 
in order to use coarse-grained methods such as DPD with soft 
potentials in a controlled way, one should use (at least) the 
persistence length as measure that should be matched between 
the coarse-grained and the atomistic models. Here, we did not 
attempt to do that systematically. 

Here, we made no attempt to determine the precise scaling 
exponents for the diffusion coefficient or the longest relax- 
ation time. Though, there are a lot of subtleties, such as the 
tube dimensions, lifetime, friction and the plateau modules, 
related to the scaling behavior 111 ill . A future publication will 
focus on them and the detailed mechanisms. 
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